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We apply recently developed techniques for simulations of moving black holes to study dynamics 
and radiation generation in the last few orbits and merger of a binary black hole system. Our analysis 
produces a consistent picture from the gravitational wave forms and dynamical black hole trajectories 
for a set of simulations with black holes beginning on circular-orbit trajectories at a variety of initial 
separations. We find profound agreement at the level of 1% among the simulations for the last orbit, 
merger and ringdown, resulting in a final black hole with spin parameter a/m = 0.69. Consequently, 
we are confident that this part of our waveform result accurately represents the predictions from 
Einstein’s General Relativity for the final burst of gravitational radiation resulting from the merger 
of an astrophysical system of equal-mass non-spinning black holes. We also find good agreement at 
a level of roughly 10% for the radiation generated in the preceding few orbits. 

PACS numbers: 04.25.Dm, 04.30.Db, 04.70.Bw, 95.30. Sf, 97.60.Lf 


I. INTRODUCTION 

Two black holes in a binary system spiral together due 
to the emission of gravitational waves. The final merger 
stage of a binary in which the black holes have compa- 
rable masses will produce a spectacular burst of gravita- 
tional radiation and is expected to be one of the bright- 
est sources in the gravitational wave sky. Mergers of 
binaries containing two stellar black holes are important 
targets for the first-generation LIGO gravitational wave 
detectors, now operating at design sensitivity in a year- 
long science data-taking run, as well as other ground- 
based detectors such as VIRGO and GEO. Knowledge 
of the waveforms from the final merger phase is impor- 
tant to improve the detectability of these sources by such 
detectors[l, 2]. Mergers of massive black hole binaries are 
important sources for the space-based LISA, currently in 
the formulation phase. Since LISA is expected to ob- 
serve these massive black hole mergers at relatively high 
signal-to- noise ratios [3], comparison of the data with cal- 
culated merger waveforms should allow a test of General 
Relativity in the dynamical, nonlinear regime. 

In the early stages of the binary inspiral, the black 
holes are widely separated and the waveforms can be cal- 
culated analytically using perturbative methods. How- 
ever, waveforms from the final merger, in which the black 
holes plunge together and form a single, highly-distorted 
black hole with a common horizon, can only be calculated 
using full 3-D numerical relativity simulations of the full 
Einstein equations. This proved be a very difficult under- 
taking and, for roughly the past decade, 3-D numerical 
relativity simulations have been beset by pernicious nu- 
merical instabililties that prevented the simulation codes 
from running long enough to evolve any significant frac- 
tion of a binary orbit. 

Recently, however, dramatic progress has been made 
in evolving the merger of equal mass binary black holes. 
Using excision to remove the singular regions within the 
horizons and a corotating coordinate system to keep the 


black holes fixed in the numerical grid, a binary has been 
evolved through a little more than an orbit, and through 
merger [4, 5], though without being able to extract gravi- 
tational waves. Simulations of excised black holes allowed 
to move through the grid on one or more orbits and then 
through merger have been carried out, with the extrac- 
tion of gravitational radiation [6]. In addition, new tech- 
niques allowing the black holes to move without the need 
for excision have been developed independently by the 
authors of this paper and another research group [7, 8]; 
these have been applied to study the final plunge, and 
merger of a binary, with the calculation of gravitational 
waves. Recently, these techniques have been applied to 
evolve a binary with nonequal masses [9] and the last 
orbit and merger of an equal mass binary [10]. It is es- 
pecially noteworthy that this progress is occurring on a 
broad front, by several independent groups using differ- 
ent techniques. 

Several major open questions in the area of binary 
black hole mergers center on the dependence of the re- 
sulting gravitational waveforms on the initial data. In 
order to use numerical relativity simulations to compare 
with data from gravitational wave detectors, we need to 
model astrophysically realistic binary black hole config- 
urations. For non-spinning equal mass black holes, how 
strongly do the gravitational waveforms for the merger 
depend on the initial data? What are the effects of spin 
and non-equal masses on the resulting waveforms? The 
answers to such questions can only be approached using 
an evolutionary analysis, in which different initial data 
sets are evolved and the resulting waveforms are com- 
pared. 

In this paper, we take a step towards answering these 
questions by evolving several initial data sets for non- 
spinning equal mass black holes through the final few or- 
bits, plunge, merger and ringdown. To accomplish this, 
we apply our new methods that allow puncture black 
holes to move freely across a grid. Using adaptive mesh 
refinement, we can resolve the dynamical regions near 
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the black holes (having length scales ~ M, where M is 
the total mass of the system and we use c = G = 1) 
and the outer regions where the gravitational waves are 
extracted (having length scales ~ (10— 100)M). Putting 
the outer boundary typically at t « 7 68M, causally 
disconnected from the dynamical regions and wavezone 
throughout most of the simulation, we can evolve stably 
for t > 800 M. The initial data sets are chosen to be 
puncture data [11]. 

Our study focuses on simulations beginning from a set 
of four inspiralling black hole configurations modelled as 
described in section II and using techniques discussed in 
section III. In section IV we calibrate the performance of 
our numerical techniques with a resolution study of sim- 
ulations from the configuration with the shortest initial 
separation. Our main results are presented in section V 
where we comparatively study the radiation waveforms 
and black hole trajectories. We find a close relationship 
between the trajectory information and the waveforms, 
providing a consistent picture of the black hole dynam- 
ics. The radiation from the final orbit, merger and ring- 
down agrees to high precision among runs from our range 
of initial configurations, with the runs from the farthest 
initial separation producing a promising waveform ap- 
proximately through the last three orbits 


II. INITIAL DATA 

We start by setting up initial data for equal mass bi- 
nary black holes represented as “punctures’’ [11]. The 
metric on the initial spacelike slice is written in the form 
gij = Sij , where i, j = 1, 2, 3, and the conformal factor 
ip = ^'bl 4- u. The static, singular part of the confor- 
mal factor takes the form ipBh = 14- J2 n =i m n/2\f - f n |, 
where the n th black hole has mass m n and is located at 
coordinate position r n . The nonsingular function u is cal- 
culated by solving the Hamiltonian constraint equation 
using AMRMG [12]. 

The punctures are initially placed on the y— axis in the 
equatorial (z = 0) plane. We need to specify the indi- 
vidual puncture mass m, coordinate position Y , and mo- 
mentum P so that the black holes are on approximately 
circular orbits and have no individual spins (they are ir- 
rotational) . To accomplish this, we adapt the quasicircu- 
lar initial data from the QC sequence given in Ref. [13]. 
These data are based in turn on the results of Cook ([14] ; 
see also [15]), who uses an effective potential method that 
minimizes the binary system energy while holding the 
orbital angular momentum fixed, to produce an approx- 
imately circular orbit. 

The parameters and physical quantities for the simu- 
lations we ran are shown in Table I. Here, Y is the initial 
coordinate position of each puncture along the y-axis, P 
is the linear momentum of an individual puncture, and m 
is the mass of the puncture. Mo is the initial total ADM 
mass of the binary, Jo is the initial total angular momen- 
tum. L is the initial proper separation of the black holes, 
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approximated by L = J ^ y/9 yy dy , where each limit Hi 
represents a point on the Schwarzschild horizon associ- 
ated with mass rm at position r*. 


Run 

±Y 

±P 

m 

Mo 

Jo 

L/Mq 

Rl 

3.257 

0.133 

0.483 

0.996 

0.868 

9.9 

R2 

3.776 

0.119 

0.488 

1.001 

0.899 

11.1 

R3 

4.251 

0.109 

0.49 

1.002 

0.928 

12.1 

R4 

4.77 

0.101 

0.492 

1.003 

0.959 

13.2 


TABLE I: Initial data parameters and physical quantities for 
the runs considered in this paper. 

We note that our runs Rl, R2, R3, and R4 correspond 
closely to the initial data for models QC-6, QC-7, QC- 
8, and QC-9, respectively, along the QC sequence. As 
shown in Ref. [13], at these relatively wide initial separa- 
tions, the parameters for these initial data sets are close 
to those derived using post-Newtonian (PN) techniques; 
see in particular Figs. 26 - 28 in Ref. [13]. 

III. SIMULATION TECHNIQUES 

We evolve the initial data with the Hahndol code[16], 
which uses a conformal (BSSN) formulation of Einstein’s 
evolution equations on a cell-centered numerical grid. 
The basic equations are the same as given in Ref. [16], 
with the exception that the evolution equation for the 
BSSN variable P has been modified for stability as sug- 
gested in [17]. In addition, to reduce high frequency noise 
associated with refinement interfaces, we add some dissi- 
pation of the form given in [18]. 

Mesh refinement and parallelization are implemented 
in our code with the PARAMESH package [19, 20]. We 
use 4 th order centered differencing for the spatial deriva- 
tives except for the advection of the shift, which is per- 
formed with 4 th -order upwinded differencing. The re- 
finement boundary interfaces are buffered with 4 th -order- 
interpolated guard cells which, at worst, may introduce 
2 nd -order errors into second derivatives. With our cur- 
rent mesh refinement implementation, our accuracy is 
limited by spatial finite differencing error from refine- 
ment interfaces. Since we do not gain much by using 
higher order time integration, we use 2 nd -order time step- 
ping via a three-step iterative Crank- Nicholson scheme. 
Even though, overall, we expect 2 nd -order convergence, 
we have found considerable advantage in using 4 th -order 
spatial differencing over 2 nd -order spatial differencing, as 
measured by the accuracy and manifest convergence of 
the Hamiltonian constraint and other quantities. 

Traditionally, puncture black holes have been evolved 
by keeping them fixed in the grid [4, 5]. This is accom- 
plished by factoring out the singular part ip bl and han- 
dling it analytically, while evolving only the regular parts 
of the metric. As explained in detail in [7], we employ in- 
stead newly developed techniques that allow the puncture 
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black holes to move freely across the numerical grid. The 
singular part is not factored out; instead the entire con- 
formal factor is evolved. Initially, the binary is set up so 
that the centers of the punctures are not located at grid 
points (as in the traditional implementation). Taking 
numerical derivatives of ^>bl causes an effective regular- 
ization of the puncture singularity through the inherent 
smoothing of finite differences. 

Our new approach allows the punctures to move freely 
by implementing a modified version of the Gamma- 
freezing shift vector. The Gamma-freezing condition gen- 
erally improves numerical stability by evolving the coor- 
dinates towards quiescence, in accord with the physical 
dynamics. Our modification to this gauge is tailored for 
moving punctures, as mentioned in [7]. Specifically we 
use 

dtP = ( 1 ) 

and 

d t B i = dtV - fidjt 1 - r ] B i (2) 

for the shift and 

d t OL = -2 aK + P'dia (3) 

for the lapse a. The lapse equation (3) is the well-known 
“H-log” slicing condition, modified with an advection 
term as in [8]. We use the initial gauge conditions (3 — 0, 
B — 0, and a = V J-2 > similar to those recommended 
in [8]. We place the punctures in the z — 0 plane and 
impose equatorial symmetry throughout. 

We use adaptive mesh refinement to produce a numer- 
ical grid having appropriate resolution in the strong-field 
dynamical regions near the black holes and in the wave 
zone. Initially, we set up the black hole binary in a nu- 
merical domain with a box-in-box refinement structure 
having an innermost refinement of hf and subsequent 
boxes of twice that resolution. We start with two boxes 
centered on each individual black hole, and then a box 
centered on the origin that encompasses both black holes. 
Subsequent boxes centered on the origin are used to give 
up to a total of 10 refinement levels. 

During the evolution the black holes move freely across 
the grid, changing the curvature in the surrounding re- 
gion; in response, the initial grid structure is changed 
adaptively. Paramesh works on logically Cartesian, or 
structured, grids and carries out the mesh refinement on 
grid blocks. If the curvature reaches a certain threshold 
(a free parameter in our code) at one point of a block, 
that block is bisected in each coordinate direction to pro- 
duce 8 child blocks, each having half the resolution of the 
parent block. If all points in all the child blocks fall below 
the threshold, those blocks get derefined. 

The box stretching from -48M to +48M in the x- 
and ^-directions and from 0 to +48M in the 0— direction 
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is fixed in place througout all runs. 1 Boxes inside this one 
can, and generally do, change adaptively as the simula- 
tion evolves; boxes outside this one maintain their orig- 
inal locations on the grid. At the initial time, there is 
a box stretching out to 24M; as the binary evolves, this 
box generally shrinks as the black holes spiral in towards 
the center. The resolution in the region between this 
box and the next (fixed) one at ±48 M is k w = 32 hf in 
all runs. We generally extract gravitational waves on a 
sphere of radius r ex = 30 M . In the early stages of the 
run, this sphere intersects the next innermost box with 
resolution h w / 2; at later times, it is completely located 
within the region having resolution h w . 

We extract gravitational waves from our simulations 
using the Weyl tensor component ^4. Our wave extrac- 
tion techniques are based on Misner’s method [21] and are 
2 nd — order accurate. They are robust and accurate even 
when the extraction radii cross mesh refinement bound- 
aries [22]. In particular, the waveforms computed at var- 
ious extraction radii r ex are preserved up to the leading 
order 1/r scaling and show no ill effects from passing 
through one or more refinement boundaries. 


IV. CALIBRATION OF SIMULATIONS 

We have performed detailed studies of the errors 
and convergence behavior of simulations run on the 
Hahndol code using fixed mesh refinement in previous 
work [7, 16, 22, 23]. The simulations presented here dif- 
fer from this earlier work in two important ways: they 
are carried out using an adaptive mesh structure that 
changes as the binaries evolve, and they are run for sig- 
nificantly longer durations. In this section, we discuss 
the calibration tests we have carried out to verify that 
the code produces robust, reliable results in these more 
challenging regimes. 

We performed these calibration tests on run Rl, which 
we ran at three different resolutions. The initial data 
parameters for Rl are given in Table I and the simulation 
parameters for these three cases are shown in Table II. 
We use a base resolution of p = 3M/32 = 0.09375M from 
which we reach the three resolutions used in the runs. 
Note that the medium resolution case, hf = p/3, has the 


Rl cases 

low 

medium 

high 

h f 

p/2 

Pi 3 

Pi 4 

outer boundary 

±768 M 

±192 M 

±768M 

T s i m 

186M 

332 M 

291M 


1 For convenience, we use M = 1 to set the scale for the compu- 
tational grid and time. Note that the actual initial total ADM 
mass for each case is Mq ss Atf , as given in Table I. 
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TABLE II: Parameters for the low, medium, and high resolu- 
tion runs of model Rl, where p — 3M/32 = 0.0937 5M. T s i m 
is the total duration of the simulation. 

outer boundary located relatively close, at ±192M; this 
was one of the earliest runs we did. Since the simulation 
runs for a total duration T S i m = 332 M, small errors from 
this outer boundary do have time to propagate in to the 
physically interesting regions of the grid. While these 
effects are generally small and have no significant impact 
on the dynamics or wave extraction, we did use a more 
distant outer boundary at ±768M in all other runs to 
eliminate this problem. This adds only a small overhead 
to the overall cost of the simulation, due to the fixed 
mesh refinement structure used in the outer regions. 

For these three runs, we scaled the criterion for refine- 
ment and derefinement to provide as close as possible 
the same grid structure at the same physical time in the 
simulations. Since it is generally not possible to match 
the grid structures exactly, pointwise convergence tests 
are of questionable value. Instead, we calculated the Ll 
norm of the Hamiltonian constraint Ch over the grid as 
a function of time. This norm was taken over all levels 
inside the box at ±48M, which includes the wave extrac- 
tion zone. Figure 1 shows the Ll norms for these three 
runs; note that the initial growth of the constraint vio- 
lation is brought to a halt after approximately 50 Af of 
evolution and then diminishes. We attribute this behav- 
ior in Ch to a gauge wave pulse that we have observed 
leaving the source region early in the simulation. The 
gauge wave has strong high frequency components and is 
thus prone to generating differencing error and reflections 
from refinement boundaries. Cr settles down somewhat 
after the gauge wave leaves the grid, as suggested by the 
plot. The curves are scaled so that, for 2 nd — order con- 
vergence, they would lie on top of each other. As Fig. 1 
shows, we get 2 nd — order convergence (or slightly better) 
for the entire course of the run. 

The gravitational waves are extracted on a sphere of 
radius r ox = 30M. The top panel of Fig. 2 shows the 
l — 2, m = 2 component of r^ 4 for the three differ- 
ent resolution runs, where r is the estimated areal ra- 
dius of the extraction sphere [22]. The waves start out 
in phase; note in particular the agreement in the initial 
pulse around t ~ 30 M. However, as the runs proceed, 
timing differences slowly accumulate over the relatively 
long orbital time scales. By t ~ 150AT, significant dif- 
ferences have accumulated in the lowest resolution case 
(solid line) compared to the other two runs; we attribute 
this to the larger inherent numerical diffusion in this low- 
est resolution run. 

A naive convergence test, in which we take differences 
directly between the waveforms at the same simulation 
times, is shown in the bottom panel of Fig. 2; these curves 
are scaled so that, for 2 nd — order convergence they would 
coincide. Comparing the top and bottom panels, it is 
clear that the differences between these curves can mainly 
be attributed to the differences in the phases of the wave- 



FIG. 1: The Ll norm of the Hamiltonian constraint violation 
is shown as a function of time for the three different resolution 
runs of Rl given in Table II. The high resolution case is shown 
with a dashed line. The medium (bold dashes) and low (solid) 
cases are scaled so that, for 2 nd — order convergence all three 
curves would lie on top of each other. The Ll norm is taken 
over all levels of the grid inside 48M, including the wave ex- 
traction region. This figure indicates satisfactory convergence 
of the Hamiltonian constraint error in our simulations. 



FIG. 2: Gravitational waveforms and a naive convergence 
test. The labels in the top panel are as in Fig. 1. The top 
panel shows the l = 2, m ~ 2 mode of ^4 for the three differ- 
ent resolution runs of Rl. The lower panel show's the differ- 
ences between these waveforms; for 2 nd - order convergence, 
the curves would lie on top of each other. Phase differences be- 
tween the waveforms account for the large differences shown. 
When the phases are shifted appropriately, the convergence 
of the waves is more manifest, as in Fig. 3 


forms. At several times, the phases are actually off by 1 r 
so that the differences between the two waves is twice as 
large as the individual waves themselves. 

A more meaningful convergence test is shown in Fig. 3. 
Here, we have shifted the waveforms in time to adjust 
for the phase differences using the techniques described 
below. The time axis has been relabeled so that the peak 
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FIG. 3: Time-shifted gravitational waveforms and a physi- 
cal convergence test. The labels in the top level are as in 
Fig. 1. In the top panel, the gravitational waveforms have 
been shifted in time so that the peak amplitude of the ra- 
diation occurs at t ~ 0. The differences between these time- 
shifted waveforms are shown in the bottom panel; these curves 
are scaled so that they would lie on top of each other for 2 nd — 
order convergence. 

of the radiation occurs at t — 0. The top panel shows the 
/ = 2, m = 2 mode of for the low (solid), medium 
(bold dashes), and high (dashes) resolution cases; note 
that the physical parts of the waveforms (t > -125 M) 
agree beautifully. We then carried out a convergence test 
using these shifted curves; the results are shown in the 
bottom panel. Again, the curves are scaled so to lie on 
top of each other for 2 nd — order convergence. Here we see 
nearly perfect convergence for the physical parts of the 
waveforms; note in particular the much smaller vertical 
scale in the lower panel of Fig. 3. 


hi 


T ex — 20 M 

Tex = 30M 

r ex = 40M 

r ex = SOM 

p! 2 

Erad 

0.0347 

0.0343 

0.0336 

0.0324 


Jrad 

0.217 

0.218 

0.218 

0.215 

p/3 

Erad 

0.0343 

0.0345 

0.0345 

0.0343 


Jrad 

0.215 

0.223 

0.226 

0.227 

p/4 

Erad 

0.0342 

0.0344 

0.0345 

0.0344 


Jrad 

0.216 

0.224 

0.227 

0.228 


TABLE III: Values of the energy E ra d (in units of Mo) and 
angular momentum Jrad (in units of Mq) carried away by 
gravitational radiation for the R1 runs calculated for different 
extraction radii and different resolutions. 

As we discussed in [7] , the gravitational waveforms can 
be used to calculate the total energy E ra d and total an- 
gular momentum J ra d carried away by the radiation. We 
calculate dE/dt and dJ/dt from time integrals of all Z = 2 
and l — 3 waveform components using Eqs. (5.1) and 
(5.2) in [13]. Integrating dE/dt gives the total energy 
loss due to the radiation, E ra d. We find that the influ- 


ence of higher modes of the waves contribute < 1% to 
the energy and even less to the angular momentum. 

The total radiated energies calculated from the wave- 
forms extracted at different radii and with different reso- 
lutions are reported in Table III. For lower resolution, the 
radiated energy depends somewhat on the extraction ra- 
dius, decreasing with increasing radius. For the medium 
and high resolution cases, however, the values are almost 
independent of the extraction radius. This indicates that 
the further out, lesser refined regions of the grid produce 
significant dissipation of the waves in the low resolution 
case, whereas when the resolution is high enough, the 
lesser refined regions do not have such an effect. These 
considerations indicate that the best radius to extract 
the energy is at r ex = 30 M, in a region refined enough 
that the energy does not significantly dissipate in the low 
resolution run, yet far enough way from the source that 
it changes only minimally for higher resolution runs. 

As shown in Table III, the radiated angular momentum 
Jrad varies by a few percent between r ex = 30M and 
r ex = 50M at high resolution. This observation seems 
to agree with the notion that the angular momentum 
depends more strongly on the longer wavelength parts 
of the waves [13], which should be extracted at greater 
distances. For all the runs reported in this paper, we use 
an extraction radius of r ex = 50 M to calculate J ra d • 


V. RESULTS 

In this section we comparatively analyze the results 
of simulations based on our R1-R4 black hole configura- 
tions. Recall that the initial data for each of these sim- 
ulations starts the black holes at different separations, 
with proper separations varying from 9. 9 Mo to 13. 2 Mo, 
and provided with sufficient angular momentum that the 
runs are estimated to be on initially circular orbits. 

We can think of each of these simulations as an approx- 
imate representation of the late-time portion of an ideal 
spacetime which begins with arbitrarily well-separated 
black holes on an inspiraling trajectory which asymptoti- 
cally approaches circular orbits. As such representations, 
there are several limitations which the initial data mod- 
els may have. In particular, the initial parameters will 
only approximate the ideal trajectory since the angular 
momentum as a function of radius may differ from the 
value required for an idealized circular orbits. Likewise, 
at finite radius, the ideal spiral trajectory can only be 
approximated by our initially circular configurations. 

Furthermore the manner by which we have mapped 
from these trajectory specifications to actual intial data 
values necessarily requires making some suppositions. 
For instance, our puncture data prescription, like almost 
all field prescriptions, contains no representation of prior 
radiation generated before the time at which the initial 
data are posed. Though it is generally expected that 
the significance of such limitations on the final merger 
simulations should be reduced if the black holes begin 
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sufficiently far apart, there is no clear way to assess just 
how significant such effects will be on the results, includ- 
ing the gravitational waveforms, before carrying out the 
evolutions. 

Our simulations, begun with varying initial separa- 
tions, should be affected by any initial modeling error 
in varying amounts, but should agree to the degree that 
they represent the ideal astrophysical spacetime. A key 
objective in our analysis is to identify universal charac- 
teristics among the different runs which, we reason, are 
then likely to correctly represent those aspects of the as- 
trophysical equal- mass binary black hole merger space- 
time. 


A. Overview of Simulations 

Our comparative analysis covers four simulations la- 
beled R1 to R4 in Table I. We evolved them all using 
the medium resolution of hf — p/3 except for Rl, where 
we used the higher hf = p/4 resolution. In all runs we 
used an initial grid setup and adaptive mesh refinement 
as described in Sec. III. We evolved all the runs to well 
after the wave signal had passed the extraction region. 
The actual amount of time is noted as T s i m in Table IV. 
For the time-slicing condition used in our simulations, 
the region where the lapse satisfies the condition a = 0.3 
corresponds roughly with the apparent horizon location. 
We thus used the moment when the two a — 0.3 regions 
around the black holes merge to specify a merger time 
Tmerge- The number of orbits for each run, iV or bits, was 
estimated from the trajectories shown in Fig. 4 and is 
taken up to the point at which the merger occurs. 



Rl R2 R3 R4 

L/Mo 

9.92 11.1 12.1 13.2 

h f 

p/4 p/3 p/3 p/3 

Tsim 

421Af 531 M 530 M 850 M 

^merger 

160Af 234 M 396 M 513 M 

-Vorbits 

1.8 2.5 3.6 4.2 


TABLE IV: Simulation parameters and results. T mer gor is the 
time at which the merger occurs, starting from the initial time 
in each run. 

A graphical overview of our four simulations is pre- 
sented in Fig. 4 showing the paths traced by the black 
hole punctures on the computational domain. These were 
obtained by numerically integrating the equation of mo- 
tion x punc = —f3{x pu nc)i which analytically gives the ex- 
act trajectory of each puncture [8]. P{x punc ) was inter- 
polated between grid points as required. 

For clarity, Fig. 4 shows only the track of one of the two 
black holes from each simulation. We have oriented each 
trajectory according to a physical reference discussed in 
Sec. VB, so that they superpose at the radiation peak, 



x / M 

FIG. 4: Paths of black holes starting from different initial 
separations. The paths are very similar for approximately 
the last orbit, indicating that the black holes follow the same 
tracks. The point of merger (estimated by a single connected 
isosurface of a = 0.3) is indicated by as an asterisk in the 
plot. 


which occurs very near the end of the puncture trajec- 
tory. R4 has the widest initial separation and completes 
the largest number of orbits. Each of the other cases, 
after an initial transient period of approximately one or- 
bit, nearly locks on to the R4 trajectory. For the final 
orbit, all four trajectories are very nearly superposed. In 
Sec. V C we study the quant iative features of these tra- 
jectories in more detail. 

Fig. 5 shows one polarization component of rT 4 for the 
runs Rl - R4, extracted at r ex = 30 M and shifted in time 
as described below. Notice that beyond about t = — 50M 
the waveforms superpose sufficiently well that it is not 
possible to distinguish the curves in the plot. The very 
strong agreement among our simulations on this part of 
the wave gives us confidence that our waveform accu- 
rately represent the astrophysical merger-ringdown signal 
to high precision. The inset shows that the agreement re- 
mains generally good going back to the beginning of each 
simulation, with the R3 and R4 runs agreeing fairly well 
over some 450Af. 


B. Gravitational Radiation 

The first step in comparing our runs is to calculate 
the energy and angular momentum carried away by the 
gravitational radiation generated by the binary system in 
our simulations. Following the discussion in Sec. IV, we 
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FIG. 5: Waveforms from runs R 1 - R4. The first 500 M are 
shown in an inset to emphasize the initial data transient part. 

measure the radiation energy extracted at r ex = 30 M, 
estimating that these will be accurate within ~ 1 %. We 
subtract the radiation energy E ra d from the initial mass 
Mq, given in Table I, to determine the final black hole 
mass in each simulation, M/ = Mo - E ra{ This pro- 
vides a physical scaling which we use to compare the 
R1-R4 simulations. Similarly, we measure the angular 
momentum content of the radiation J ra d as extracted at 
r ex = 50M, estimating the accuracy to be within a few 
percent. Subtracting this from the initial ADM angular 
momentum Jo, we calculate the spin parameter of the 
final black hole a/Mf = (Jo - J rac i)/M^. The results are 
summarized in Table V. 



Er ad /M J 

ra d/M 2 

a/M 

Mqn /M agAr/M 

Rl 

0.0356 

0.246 

0.694 

1.005 

0.721 

R2 

0.0369 

0.272 

0.691 

1.002 

0.686 

R3 

0.0381 

0.306 

0.689 

1.004 

0.694 

R4 

0.0387 

0.325 

0.702 

1.004 

0.693 


TABLE V : Energy and angular momenta for the radiation and 
final black hole. E ra d and J ra d are measured at r ex = 30M, 
and r ex = 50M, respectively. Mqn and clqn are calculated 
independently from the quasi-normal fits of the ringdown 
waveforms, and agree well with the values deduced from the 
radiative losses. 

We note that energy and angular momentum content 
of the radiation is almost entirely contained in the l = 2 , 
m = ±2 spin - 2 -weighted spherical harmonic compo- 
nents, with other components entering at the 1 % level. 
In the remainder of our waveform analysis we concentrate 
exclusively on the leading component, 

rip4(0,ip) ~ r^4(22) (-2^22(0, 0) + -2^2 -2 (#> <£)) (4) 

where, for simplicity, we have suppressed retarded time 
dependence. Hereafter we suppress the multipole labels 
and refer to the leading component simply as r$ 4 . The 



FIG. 6 : Amplitudes (absolute value of (complex) ^ 4 ) of the 
waves. The curves have been shifted such that the maxima 
are all at time 0. The inset zooms into the peak showing the 
strong agreement from t = —50 M on. 


two polarization components of the radiation are rep- 
resented in the real and imaginary parts of 7 *^ 4 . We 
follow Refs. [13, 24], representing our waveforms by 
r'F 4 = Aexp(— i<p po i), where the amplitude A and po- 
larization phase ippoi, like r4> 4 , are functions of time. Of 
course, any complex time-series could be represented in 
this form, but it is particularly valuable if the radiation 
exhibits a circular polarization pattern, meaning that A 
and (ppoi vary slowly compared to the wave frequency 
timescale. Such radiation will be circularly polarized to 
an observer on the system’s rotational axis, varying to 
linearly polarization for an observer on the equatorial 
plane. We find that the radiation produced in our simu- 
lations shows strong circular polarization, which we will 
utilize in comparing the radiation from our four simula- 
tions. 

As we are interested in the degree to which our vari- 
ous runs may be taken to be different models of the same 
physical merger spacetime, we must define a physical ba- 
sis for comparing them. As physically expected, runs 
beginning with more separated black holes take a longer 
time to reach the point of merger. For each run, a phys- 
ical reference time is recognizable by the point at which 
the radiation reaches its peak amplitude; we define t = 0 
at this point for our comparisons. 

Fig. 6 shows the wave amplitudes A(t) from all four 
runs. Through the strong radiation peak after t = —50 M 
all four wave amplitudes show striking universality in 
the compared results, with agreement among all runs 
to about 1%. This period of strong agreement cov- 
ers roughly the last orbit, plunge and ringdown of the 
merger. Agreement within about 10% is maintained 
among the R2-R4 simulations for most of each run with 
slightly more difference in the R 1 run. The smooth shape 
of the peak, lacking any sign of the several wave cycles 
spanning the peak (see Fig. 5) is an indication of the 
circular polarization pattern discuss above. The only 
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FIG. 7: The phase angle vs. time. The phase is made to 
agree at time OM. But after that the phases agree very well 
showing that the ringdown is the same for all the runs. And 
it agrees before that time showing that the black holes follow 
the same dynamics independent of their initial separation. 


clearly non- universal feature in the wave amplitudes is 
a small and brief burst lasting about TOM at the begin- 
ning of each run. We interpret this burst as “spurious 
radiation” content in the initial data. As generally ex- 
pected the amplitude of the spurious radiation lessens as 
the initial separation of the black holes is increased. 

Since the black holes in our various runs (which all be- 
gin on the y- axis) undertake differing amounts of orbital 
motion before merger, we must expect differing orienta- 
tions for the systems at the point of merger. For mean- 
ingful comparsions we must rotate the data from each 
system to align them with repect to some physical ref- 
erence. As our reference, we will orient the systems so 
that the polarization phase <p po i of the radiation passes 
though zero at the moment of peak amplitude t = 0. 
We use this orientation for all figures with wave phase 
information, and in the trajectory comparison in Fig. 4. 

The polarization phase (p po i is shown in Fig. 7. At 
time OMo the phases are set to agree. However, they 
keep agreeing after that, showing that the ringdown fre- 
quency is the same. The plot shows generally good phase 
agreement among all runs, to within a small fraction of 
a wave cycle, except for a brief period at the beginning 
of each run when the radiation is dominated by spuri- 
ous radiation associated with initial data modeling error, 
and is not circularly polarized. For circularly polarized 
radiation it is meaningful to define an instantaneous fre- 
quency for the wave, u = d<p po i/dt , shown in Fig. 8, 
which allows a more detailed comparison of the simula- 
tion waveforms. Clearly can be seen that apart from the 
noise due to the smallness of the waves initially and at 
the end, all runs have the same frequency evolution from 
about 60Mo before the merger. Before that the frequency 
evolution compares similarly among the simulations with 
some wavy ness which may correspond to some ellipticity 
in early part of each simulation’s insprial. 



FIG. 8: Waveform frequency as a function of time. The agree- 
ment from t = —50 M on is astounding. 

The final (constant) frequency is the frequency of the 
ringdown. Fitting an exponential decay to the ampli- 
tudes and using this final frequency we estimate, using 
the procedure in Ref. [25], the mass and angular mo- 
mentum of the black hole formed in the merger. These 
estimates provide a description of the final black hole as 
determined by its perturbative dynamics. The results 
are listed in Table. V where they can be compared with 
the independent estimates obtained by subtracting the 
radiation losses from the initial values. The good agree- 
ment, to one or two percent, provides another measure of 
energy and angular momentum conservation in each sim- 
ulation. We also note the strong agreement among our 
simulations from different separations on the final state 
of the remnant black hole, with every measure consistent 
with the same value, a/M — 0.69(±2%). 

C. Trajectory analysis 

In this section we consider the particle- like dynamics 
of our simulations defined by the coordinate trajecto- 
ries of our black hole punctures. It is important to be 
careful in interpreting such coordinate-dependent infor- 
mation which may include non-physical gauge features. 
Nonetheless it is worth noting that the T-freezing gauge 
condition applied in our simulations, to the extent that it 
approximates T l — Jj =0 is similar to the Dirac gauge 
applied in some post-Newtonian calculations, and might 
be expected to provide a sensible coordinate system at 
least in weak field regions. In any case, since we apply 
similar coordinate conditions in each run, the coordinate- 
based puncture trajectories provide another opportunity 
to identify universal features in our simulations from dif- 
ferent initial separations. 

Refer again to Fig. 4, which shows the tracks of one 
of the punctures from each of the runs R1 - R4 as they 
spiral into the center. As the black holes descend deeper 
into the strong field, we see that their paths lock on to 
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FIG. 9: The coordinate separation between the punctures is 
shown for runs R1 * R4 as a function of time. 


a universal trajectory that takes them into the plunge 
and subsequent merger. This behavior can also be seen 
in Fig. 9, which shows the coordinate separation between 
the punctures, r = ( x 2 + y 2 ) 1 / 2 as a function of time. In 
this case, there is strong agreement among the runs after 
t & -50 M. In the earlier part of the simulations there 
are clear differences among the runs, which are suggestive 
of ellipticity in the initial orbital motion. 

One way to estimate the utility of such coordinate in- 
formation is by comparison with our much more invariant 
waveform data. We can, for instance, compare the coor- 
dinate orbital frequency Cl with the waveform frequency 
examined above. With suitable coordinates, in the weak- 
field limit, we would expect the orbital frequency to be 
approximately equal to half the gravitational wave fre- 
quency. In Fig. 10 we compare the coordinate orbital fre- 
quency Cl with half the gravitational wave frequency from 
our R4 run. We find, that if we shift the orbital freqency 
data by about 33 M 0 , the two curves match very well. 
Despite some noise in our wave frequency, most features 
in the orbital frequency are tracked in the wave frequency 
as well. Exceptions occur at the very beginning of the 
simulation, where our coordinates necessarily start with 
the punctures non-moving (hence Cl = 0) before the brief 
period through which the coordinate puncture velocity 
seems to adapt well to the physical dynamics. Likewise, 
at late times, we see that the orbital puncture frequency 
continues to grow while the radiation frequency saturates 
at the quasinormal ringing frequency. This seems to cor- 
respond to the expectation that, at late times, the coordi- 
nate motion of the punctures decouples from the process 
of radiation generation, with the punctures continuing to 
fall into the newly formed black hole while gravitational 
ring-down radiation is generated in final black hole’s per- 
tubative potential barrier region somewhat outside the 
horizon. 

The 33Mo time shift required to realize this agreement 
can be interpreted as the time required for the gravi- 
tational radiation to propagate out to point where it is 



FIG. 10: Orbital and radiation frequencies for R4. We show 
the evolution of the orbital frequency Cl, calculated from the 
coordinate motion of the punctures, compared with the wave 
frequency. The good corresondence motivates a closer exam- 
ination of the coordinate-based puncture trajectories. (Note 
that the radiation frequency is divided by 2 since uj = 2fh) 



FIG. 11: The frequencies over times. Shown are frequencies 
calculated from the waves showing the agreement of the differ- 
ent runs and two waveforms calculated using Post-Newtonian 
approximations. 

extracted at R = 30 Mq. We can utilize this time corre- 
spondence to roughly associate phenomena occurring in 
the strong field region of the simulations with features 
in the radiation. We have used this time- shift to com- 
pare the time at which peak radiation is generated in our 
simulations with T me r 5 e in Table IV. In our simulations, 
we find that the merger occurs about 12M before the 
radiation peak is generated. 

In Fig. 11 we compare the orbital frequencies of the 
puncture trajectories among our R1-R4 runs. As was the 
case with the wave frequencies we see excellent agreement 
among the runs after t = -50 M. Here we have focused 
on the earlier region where the results are not quite as 
univeral. In each run we notice that after a period of ini- 
tial angular acceleration the orbital frequencies agree to 
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within a few percent. This initial transient period is at 
least partially a coordinate effect (since Cl must begin at 
0), but the strong similarity with wave-frequencies sug- 
gests that there may be a physical basis for the descrep- 
ancies as well. Particularly for R3 and R4, the shape of 
the curves suggests a slow oscillation, perhaps some ellip- 
ticity in the motion. For an external comparison, we have 
included the post-Newtonian frequency evolution 1.5PN 
and 2PN order, as provided in Ref. [26], positioned in 
time so that all curves agree near t = —200 M. The cor- 
respondence with these PN curves is certainly striking, 
though we caution that the 2.5PN curve does not agree 
so well, with its peak frequency topping out too low for a 
useful comparison here. In subsequent work we intend to 
explore comparisons with post-Newtonian calculations in 
more detail. 


VI. DISCUSSION 

We have presented a set of numerical simulations rep- 
resenting the last few orbits and merger of an equal-mass 
non-spinning binary black hole system. Though the ini- 
tial data differ, with the black holes starting out at sep- 
arations ranging from 9M to 13M, the calculated wave- 
forms are dominated by universal characteristics. Over 
the period beginning from about 50 M before the gravi- 
tational wave peak, and covering approximately the last 
orbit, all our simulations show profound agreement, with 
differences among the waveforms at no more than about 
1%. Consequently, we are confident that this part of 
our waveforms accuartely represents the predictions from 
Einstein’s General Relativity for the final burst of gravi- 
tational radiation resulting from the merger of an astro- 
physical system of equal-mass non-spinning black holes. 

We also see good agreement among our simulations in 
the gravitational radiation generated in the several pre- 
ceding orbits which we have simulated. Excepting an ini- 
tial transient period about 100M for each run, the earlier 
portion of the waveforms show good agreement, within 
approximately 10% in phasing and amplitude. For our 
longest simulation (R4) this suggests that we have pro- 
duced a good approximation of the astrophysical wave- 
form prediction covering more than 400M before the ra- 
diation peak. We are further encouraged by good agree- 
ment over much of this period among the frequency evo- 
lution of our simulated waveforms (especially R4), with 
the second-order (2PN) post-Newtonian predictions. We 
are currently exploring the correspondence of our wave- 
forms with post-Newtonian calculations, to be developed 
in more detail in a future publication. 

That the late- time part of the simulations (the final or- 
bit and thereafter) shows stronger universality than the 
earlier part of the simulations supports idea that the late 
dynamics of the system is dominated by the strong inter- 
action of the holes, and radiative losses, which have the 
effect of reducing dependence on prior conditions. For 
the remnant black hole formed in the merger, our sim- 


ulations consistently predict a spin within about 1% of 
a/m — 0.69. 

Our simulations have employed newly developed nu- 
merical relativity techniques for evolving black holes, 
which allows the black holes to propagate accurately 
across the numerical domain[7, 8]. This approach does 
not require excision of the black hole interior from the 
computational domain. We have calibrated our approach 
on one of our black hole configurations, R1 , demonstrat- 
ing 2 nd -order convergence, and waveform accuracy at the 
1% level in the last orbit. For all our simulations we have 
found good energy and angular momentum conservation 
as measured by comparing the mass and spin of the final 
black hole, as measured by its quasinormal ringing with 
the expected remainder after radiative losses. In future 
work we expect to continue to refine our techniques for 
accuracy and efficiency over long-lasting simulations. 

We have also studied the trajectories traced out by the 
motion of the black holes in our numerical coordinate sys- 
tem. These trajectories suggest a coherent picture of the 
system’s evolution which is qualitatively, and in some 
ways, quantitatively consistent with invariant informa- 
tion measured in the radiation waveforms. This corre- 
spondence provides some foundation for giving a tenta- 
tive physical interpretation to some coordinate-based in- 
formation in our simulations, such as the number and 
rate of inspiral orbits and recommends further research 
toward a deeper understanding the properties of the co- 
ordinate systems which we have applied in these simula- 
tions. 

Our comparative analysis provides some insight into 
the quality of the initial data models which we have ap- 
plied. We have seen an indication of a small amount 
(with negligible energy content) of spurious radiation 
arising from initial modeling error. As generally expected 
the scale of this spurious radiation decreases for increas- 
ingly well-separated initial data. For sufficiently long- 
lasting runs this spurious radiation is well-segregated in 
time, limiting its direct significance in interpreting the 
merger radiation. Of perhaps greater concern are indi- 
cations which suggest ellipticity in the inspiral trajec- 
tories and waveforms. Two aspects of our initial data 
model may be contributing to this effect. First, the will 
be some level of mismatch between the initial separa- 
tion of the black holes and the specified initial value of 
angular momentum. We comment, in this regard, that 
other simulations we have looked at in preparation for the 
work presented here hint that the early transient part of 
each simulation, as well as measures such as the merger 
times, may depend sensitively on small (1%) changes in 
the initial angular momentum. Secondly, as is custom- 
ary, the data we have applied have the black holes set on 
initially circular trajectories, with vanishing radial mo- 
mentum components. This may also be expected to lead 
to ellipticity, as has been noted in binary neutron star 
simulations [27]. We expect to explore these and similar 
concerns in future simulations. 

Taken together with other recent progress in numeri- 
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cal relativity, these results herald a new age of numer- 
ical simulations applied to further characterize and un- 
derstand strong-field binary black hole interactions and 
merger radiation. Future applications will begin to ex- 
plore the physical parameter space of these systems, to 
study, particularly, the effect on the radiation waveforms 
of individual black hole spin, and variations in the black 
hole mass ratio. 
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